# ===============================================
# Paris ratchet paper
# Code to make figure temperature summaries across treaties
# ===============================================

#### Table of contents ####

#' Outline
#' -- Load in global temperature data
#' -- Load RCP concentrations and temperatures to get the linear CO2-temperature relationship
#' -- -- rcp_85; co2; rcp_co2
#' -- Then regress temperature on concentrations, and save alpha and beta
#' -- Then paste in scenario summaries and get consistent 90% CIs
#' -- Then make calculations for Kyoto temperatures and CIs
#' -- And plot


#### Initialize ####

library(tidyverse)
library(readxl)
library(broom)

## Set working directory


#### Get global annual temperatures

## NOT RUN
# read_csv("https://www.ncei.noaa.gov/access/monitoring/climate-at-a-glance/global/time-series/globe/land_ocean/12/1/1850-2023/data.csv",
#                     skip = 4) %>%
# select(year = 1, temp_anomaly = 2) 

temp_noaa <- read_csv("input/temp_noaa.csv") %>% 
  ## Rebase temperature
  mutate(period_1850_1900 = ifelse(year %in% seq(1850, 1900, 1), "pre", "post")) %>% 
  #' NOAA uses 1901-2000 baseline, so re-base to 1850-1900
  group_by(period_1850_1900) %>% 
  mutate(mean_period_temp = mean(temp_anomaly)) %>% 
  ungroup() %>% 
  mutate(mean_period_temp = replace(mean_period_temp, period_1850_1900 == "post", NA),
         pre_period_temp = max(mean_period_temp, na.rm = T),
         temperature = temp_anomaly + pre_period_temp) %>% 
  select(year, temperature)


#### Get RCP8.5 data to have some temperatures
rcp_temp <- read_excel("input/rcp_scenarios.xlsx", sheet = "rcp_temperature_projections") %>% 
  select(year, rcp_85_50) %>% # Here we have median and 1sd of RCP8.5
  # Re-base the -0.61 from 1850-1900
  mutate(rcp_85_50 = rcp_85_50 + 0.61) %>%  #; note 0.61 correction comes from this sheet
  # Get values of temperature and CO2 in 2010 (observed) and 2100 (projected)
  filter(year %in% c(2010, 2090)) %>% # 2010 and 2090
  select(rcp_85_temp = rcp_85_50)


#### Get CO2 concentrations
rcp_co2 <- read_excel("input/rcp_scenarios.xlsx", sheet = "rcp_co2_concentrations") %>% 
  filter(year %in% c(2010, 2090:2100)) %>% 
  mutate(time_period = ifelse(year == 2010, "Observed", "Projected")) %>% 
  group_by(time_period) %>% 
  summarize(rcp_85_concentration = mean(rcp_85_concentration)) %>% 
  mutate(year = c(2010, 2100))


#### Linear relationship between CO2 and temperature
#' https://www.nature.com/articles/nature08047
true_alpha <- bind_cols(rcp_temp, rcp_co2) %>% 
  lm(rcp_85_temp ~ rcp_85_concentration, data = .) %>% 
  tidy() %>% 
  filter(term == "(Intercept)") %>% 
  pull(estimate)

true_beta <- bind_cols(rcp_temp, rcp_co2) %>% 
  lm(rcp_85_temp ~ rcp_85_concentration, data = .) %>% 
  tidy() %>% 
  filter(term == "rcp_85_concentration") %>% 
  pull(estimate)


#### Paste in scenario summaries data
scenarios <- data.frame(
  stringsAsFactors = FALSE,
            author = c("rogelj2016", "rogelj2016", "meinshausen2022"),
            source = c("table1", "table1", "abstract"),
          scenario = c("No policy", "Paris", "Glasgow"),
       temperature = c(4.1, 2.7, 1.95),
          lower_ci = c(3.1, 2.1, 1.4),
          upper_ci = c(4.8, 3.2, 2.8),
          range_ci = c("80% CI", "80% CI", "90% CI"),
           range_z = c(1.28, 1.28, 1.645)) %>% 
  mutate(lower_se = abs((temperature - lower_ci)/range_z),
         upper_se = abs((upper_ci - temperature)/range_z),
         lower_90 = temperature - 1.645*lower_se,
         upper_90 = temperature + 1.645*upper_se) %>% 
  mutate(co2_concentration = (temperature - true_alpha)/true_beta)
scenarios


#### Now get Kyoto values
#' Wigley thinks KP saves 50ppm CO2
kyoto_temp <- true_alpha + true_beta*(scenarios$co2_concentration[scenarios$scenario=="No policy"]-50)
kyoto_temp_diff <- scenarios$temperature[scenarios$scenario=="No policy"] - kyoto_temp
kyoto_lower_90 <- scenarios$lower_90[scenarios$scenario=="No policy"] - kyoto_temp_diff
kyoto_upper_90 <- scenarios$upper_90[scenarios$scenario=="No policy"] - kyoto_temp_diff


#### Combine Kyoto calculations with those in the summaries, and get 90% CI  
plot_scenarios <- scenarios %>% 
  add_row(author = "wigley", source = "sr_calc", scenario = "Kyoto",
          temperature = kyoto_temp, 
          lower_90 = kyoto_lower_90,
          upper_90 = kyoto_upper_90) %>% 
  select(scenario, temperature, lower_90, upper_90) %>% 
  mutate(year = case_when(scenario == "No policy" ~ 1992,
                          scenario == "Kyoto" ~ 1997,
                          scenario == "Paris" ~ 2015,
                          scenario == "Glasgow" ~ 2021)) 


#### Define theme
gg_theme <- theme(legend.position = "none",
                  panel.background = element_rect(fill=NA), 
                  panel.grid.major.x = NULL,
                  panel.grid.major.y = NULL,
                  panel.ontop = FALSE,
                  text=element_text(size=12),
                  axis.title = element_text(family = "serif"),
                  axis.text = element_text(family= "serif"),
                  axis.line = element_line(size = 1, colour = "black"))

#### Make plot
temperature_plot <- plot_scenarios %>% 
  ggplot(., aes(year, temperature)) +
  coord_cartesian(xlim = c(1990, 2030))+
  scale_y_continuous(limits = c(0, 5.5),
                     breaks = c(0, 0.5, 1, 1.5, 2, 3, 4, 5)) +
  geom_rect(aes(xmin=1980, xmax=2100, 
                ymin=lower_90[scenario=="No policy"], 
                ymax=upper_90[scenario=="No policy"]), 
            fill = "#FFD492", alpha = 1) +
  annotate("text", y = 4.5, x = 2030, hjust = 1, 
           label = "Projected warming in 2100 without climate policy",
           family = "serif") +
  geom_rect(aes(xmin=2015, xmax=2100, ymin=1.5, ymax=2), fill = "#BFE7FF") +
  annotate("text", y = 1.75, x = 2030, label = "Paris goal", hjust = 1,
           family = "serif") +
  geom_errorbar(data = plot_scenarios %>% filter(year>1995),
                aes(x = year, ymin = lower_90, ymax = upper_90), width = 1) +
  geom_text(data = plot_scenarios %>% filter(year>1995),
            aes(x = year, y = upper_90, label = scenario,
                family = "serif"), vjust = -0.5) +
  geom_line(data = plot_scenarios %>% filter(year>1995),
            aes(year, temperature), linetype = 2) +
  labs(x = NULL, y = "Temperature increase from 1850-1900") +
  # geom_line(data = berkeley_temp %>% mutate(scenario = "1"),
  #            aes(x = year, y = temperature)) +
  geom_line(data = temp_noaa %>% mutate(scenario = "1"),
            aes(x = year, y = temperature),
            color = "black",
            size = 0.75) +
  gg_theme
  
ggsave(plot = temperature_plot, 
       filename = 'plots/ratchet_temperatures.pdf', units = 'in', height = 4, width = 7)


   
